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Perturbations of Schwarzschild-Droste black holes in the Regge- Wheeler gauge benefit from the 
availability of a wave equation and from the gauge invariance of the wave function, but lack smooth- 
ness. Nevertheless, the even perturbations belong to the C° continuity class, if the wave function 
and its derivatives satisfy specific conditions on the discontinuities, known as jump conditions, at the 
particle position. These conditions suggest a new way for dealing with finite element integration in 
the time domain. The forward time value in the upper node of the (t,r*) grid cell is obtained by the 
linear combination of the three preceding node values and of analytic expressions based on the jump 
conditions. The numerical integration does not deal directly with the source term, the associated 
singularities and the potential. This amounts to an indirect integration of the wave equation. The 
known wave forms at infinity are recovered and the wave function at the particle position is shown. 
In this series of papers, the radial trajectory is dealt with first, being this method of integration 
applicable to generic orbits of EMRI (Extreme Mass Ratio Inspiral) . 

PACS numbers: 04.25.Nx, 04.30.Db, 04.30.Nk, 04.70.Bw, 95.30.Sf 
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I. INTRODUCTION 



It is currently affirmed that the core of most galaxies 
host supermassive black holes, on which stars and com- 
pact objects in the neighbourhood inspiral and plunge- 
in. The EMRI (Extreme Mass Ratio Inspiral) sources are 
characterised by a huge number of parameters that, when 
spanned over a large period, produce a huge number of 
templates for the LISA (Laser Interferometer Space An- 
tenna) project Thus, matched filtering (the superpo- 
sition of the expected and received signals) [2| is aided by 
stochastic methods @; in alternative, other more robust 
methods based on covariance or on time and frequency 
analysis are investigated 0, 341] ■ Furthermore, if the sig- 
nal from a capture is not individually detectable, it still 
may contribute to the statistical background [t| [l(J ■ De- 
tection strategies are to be tested by the LISA simulators 

Challenges by EMRIs arise increasing interest amongst 
astrophysicists and those studying data analysis, but also 
for theorists. Relativists, especially, have been drawn 
to study EMRIs for the determination of the motion of 
the captured mass and because of the impact that back- 
action has on the gravitational wave forms (see [l4| where 



a comprehensive introduction to these topics is given). 

The captured compact object may be compared to a 
small mass m perturbing the background spacetime cur- 
vature of a large mass M and generating gravitational ra- 
diation. Perturbation methods from the seminal work of 
Regge and Wheeler, henceforth RW [HI], were an obvious 
aid for studying the case we are dealing with. Approaches 
based on numerical relativity (NR) appear not be best 
suited, nor do the post-Newtonian (pN) expansions: the 
former for the two-scale problem (one small scale for the 
neighbourhood of the particle, one large scale for the ra- 
diation at infinity); the latter for the small field and low 
velocity constraints. To EMRIs, nevertheless, Damour, 
Nagar and co-workers [l6l - [l9l | , and Yunes et at. have 
applied the EOB (effective one body) method, though 
based on the pN approach. Furthermore, both pN and 
perturbation methods are found to be in agreement in 
their common domain of applicability by Blanchet et al. 
[2TI |22| , as well as EOB and perturbation methods have 
shown their synergy thanks to Barack et al. [23| . Fi- 
nally, it is to be mentioned that NR analysis is slowly 
progressing towards unequal mass ratios [24j . 



Perturbation methods 



"This arXiv version differs from the one published in Phys. Rev. D 
for the use of British English and other minor editorial differences. 



EMRIs have revived the interest for the perturbative 
relativistic two-body problem, which was investigated 
with a semi-relativistic approach some 40 years ago by 
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Ruffini and Wheeler [25|, |26|, and then more fully by 
Zerilli [27M30| . These authors treated the source of the 
perturbations in the form of a radially falling particle. 
Zerilli's results opened the way for others to study to 
first order the orbital motion of the captured mass in a 
fully, although linearised, relativistic regime. Work done 
up until the end of the '90s involved the captured star 
moving on the geodesic of the background field and being 
unaffected by its own mass and the emitted radiation. 

We shall refer to the first solution of the Einstein equa- 
tion as the work of Droste [3lT[33l | and Schwarzschild |IH , 
collectively as SD, this double authorship having been 
justified by the historians, as recalled by Rothman f35j . 
Much research has been devoted to the study of SD black 
hole perturbations in the RW gauge, first in a vacuum 
11511 . and then by Zerilli in the presence of a source 
130( 1 . An outflow of pub lications by Davis, Press, Price, 
Ruffini and Tiomno [36l - l41| , but also by individual schol- 
ars like Chung [Hj], Dymnikova |43| appeared in the '70s. 
It was followed by a slowdown in the '80s and '90s, with 
the exception of the forerunners of the Japanese school 
as Tashiro and Ezawa (44|, Nakamura with Oohara and 
Koijma [45| or with Shibata [46|. The aim of all these 
authors was to analyse especially the amplitude and spec- 
trum of the radiation emitted by the falling mass. Sim- 
ulations were performed in frequency domain with the 
particle falling from infinity. The numerical accuracy was 
recently improved by Mitsou (4?J . It was in 1997, that a 
fall from finite distance was analysed by Lousto and Price 
[48| . The latter two authors were the first to use a finite 
difference scheme in time domain [49| . Following on this 
work, Martel and Poisson [Ho| were able to parametrise 
the initial conditions reflecting past motion of the parti- 
cle and resulting into an initial amount of gravitational 
wave energy. All of this work was carried out in the RW 
gauge. 

It is less than 15 years since we began to have methods 
for evaluating the back-action for point masses in strong 
fields thanks to two concurring situations. On one hand, 
the theorists progressed in understanding radiation reac- 
tion and obtained formal prescriptions for its determina- 
tion. The solution was brought by the self-force equation, 
first by Mino, Sasaki and Tanaka [Hl[, Quinn and Wald 
[5^ , and later by Detweiler and Whiting [53| , all various 
approaches yielding the same formal expression. On the 
other hand, for the above mentioned detection of cap- 
tures, the requirements from LISA posed constraints on 
the tolerable amount of phase-shift on the wave forms, 
caused by radiation reaction. This impulse - project ori- 
ented - compelled the researchers to turn their efforts in 
finding an efficient and clear implementation of the pre- 
scription made by the theorists. The development was 
again pursued in perturbation theory. This time, though, 
the small mass m is correcting the geodesic equation of 
motion on a fixed background via a factor 0(m). 

Actually, the advances in the field by e.g. Warburton 
and Barack [54[ and in radiation gauge by Kcidl et al. l55ll 
have recently tackled the self- force in Kerr geometry [56( . 



Such initial results in Kerr geometry do not mean that 
all major issues in non-rotating black holes have been 
solved. For instance, the self-consistent evolution of an 
orbit, along the lines suggested by Gralla and Wald (57j . 
is not yet applied even to the simpler case of non-rotating 
black holes. 



B. A different integration approach 

The complexity in assessing the continuity of the per- 
turbations (and the associated numerical computation of 
derivatives) at the position of the particle, has led Barack 
and Lousto [58j, among other motivations such as the 
compatibility of the self-force to the (Lorenz-FitzGcrald- 
de Donder) harmonic gauge [59l - l6l1 |. to convey their ef- 
forts to this gauge. This alternative enterprise is by no 
means priceless, the price being the unavailability of a 
wave equation like in the RW gauge and of the gauge 
invariance of the wave function, as determined by Mon- 
crief (62). In this context, the recent investigations of a 
gauge independent form of the self-force (63j may change 
sensibly the trade-off on the gauge choice. 

The choice of a gauge should be suited to the type 
of problem under scrutiny and thus we do not adopt a 
position on which gauge should be preferred a priori. 
We find of some interest to revert a weakness of the RW 
gauge, namely the lack of smoothness, into a building 
tool of a new integration method. 

Indeed, we propose herein a finite element method of 
integration based on the jump conditions that the wave 
function and its derivatives have to satisfy for the SD 
black hole perturbations to be continuous at the position 
of the particle. We first deal with the radial trajectory 
and the associated even parity perturbations, while in a 
forthcoming paper [64[ we shall present the circular and 
eccentric orbital cases, referring thus to both odd and 
even parity perturbations. 

Herein, we present a first order method and show that 
it suffices to acquire a well-behaved wave function at in- 
finity and at the particle position. However, a higher 
order scheme has been derived and tested [65| . 

The main feature of the indirect (-I-) method consists 
in avoiding the explicit integration of the wave equation 
(the source term with the associated singularities and the 
potential) whenever the grid cells are crossed by the par- 
ticle. Indeed, the information on the wave equation is 
implicitly given by the jump conditions. Conversely, for 
cells not crossed by the particle, we retain the classic ap- 
proach for integrating wave equations [49|, Ho| , henceforth 
named LPMP. 

For the computation of the back-action, the -I- method 
ensures a well behaved wave function at the particle po- 
sition, since the approach is governed by the theoretical 
values of the jump conditions. 
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C. Structure of the paper 

The structure of the paper is as follows. Section II is 
a brief reminder of the formalism on perturbations for a 
SD black hole in the RW gauge. In Sec. Ill, we give an 
overview about the jump conditions on the wave function 
and its derivatives. In Sec. IV, we introduce the scheme 
by detailing one of the four physical cases that a particle 
crossing the mesh cell during its infall may encounter, 
while for the other three similar cases we simply provide 
the final result. In Sec. V some wave forms at infinity are 
shown, discussed and compared to those obtained with 
the LPMP method. Furthermore, the theoretical and 
numerical jump conditions at the position of the particle 
are displayed and commented. In the conclusions, Sec. 
VI, the perspectives are addressed. 

Geometric units (G = c = 1) are used, unless stated 
otherwise. The metric signature is (—,+,+,+). 

II. BLACK HOLE PERTURBATIONS 

The Zerilli-Moncrief (ZM) equation rules even-parity 
waves in the presence of a source: a freely falling point 
particle, that generates a perturbation for which the dif- 
ference from the SD geometry is small. The energy- 
momentum tensor T^ v is given by the integral of the 
world line of the particle, the integrand containing a four- 
dimcnsional invariant 8 Dirac distribution for the repre- 
sentation of the point particle trajectory. The vanishing 
of the covariant divergence of T^ v is guaranteed by the 
world line being a geodesic in the background SD geom- 



etry, represented by the metric tensor g^. Finally, the 
complete description of the emitted gravitational waves 
is given by the symmetric tensor h^ v <C g^v 

The formalism can be summarised as follows. Any 
symmetric covariant tensor can be expanded in spherical 
harmonics [29[ . Because of the spherical symmetry of the 
SD field, the linearised field equations are cast in the form 
of a rotationally invariant operator on h^ u . This term is 
equated to the energy-momentum tensor, also expressed 
in spherical tensorial harmonics. That is: 



oc r M1/ [d(r„)] , (1) 

where the S(r u ) Dirac distribution represents the point 
particle unperturbed trajectory r u . The rotational in- 
variance is used to separate out the angular variables in 
the field equations. For the spherical symmetry of the 
2-dimcnsional manifold on which t, r arc constants under 
rotation in the 6,<f> sphere, the ten components of the 
perturbing symmetric tensor transform like three scalars, 
two vectors and one tensor: 




In the Regge- Wheeler- Zerilli formalism, the source 
term for the odd perturbations vanishes for the radial 
trajectory and, given the rotational invariance through 
the azimuthal angle, only the index referring to the po- 
lar or latitude angle survives. The even perturbations, 
going as (—1)', are expressed by the following matrix: 



huu 



2f ) H l Y l 


H[Y l 


h' Y> e 


sym 




h[Y< e 


sym 


sym 


r 2 [K'Y 1 +G l Y* gg ] 


sym 


sym 


sym r 2 sin 2 



r 2 G l 
K l - 



(YL- cot 9Y 



<:-) 



(2) 



G 



Y 



+ cot 0Y l e 



where H l , H[ , H l 2 , h l , h\ ,K l ,G l arc functions of (t, r) and 
the spherical harmonics Y l obviously of (9,<fr). Two for- 
malisms describe the evolution of the perturbations in 
terms of a single wave function and a single wave equa- 
tion. One, due to Moncrief (6^1, gauge invariant and 
developed about a 3-geometry of a t = constant hyper- 
surface, refers solely to the Hi, K l ,G l ,h[, perturbations. 
The other, due to Zerilli (234131, uses the RW g au § e 



G 



h[ = 0. In the RW gauge, the Moncrief 



wave function is given by: 



A + I 



Xr + 3M 1 2 



pK 
dr 



i \ 1 



, (3) 



where the Zerilli |28[ normalisation is used for The 
dimension of the wave function is such that the energy is 
proportional to J °° \E' 2 dt. Finally, the wave equation is 
given by: 



dr 



> 2 



dt 2 



V l (r)^(t,r) = S t (t,r) , (4) 



where r* = r + 2M\n(r/2M — 1) is the tortoise coordi- 
nate. The potential V (r) is expressed by: 
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, / 2M\ 2A 2 (A+l)r 3 +6A 2 Mr 2 + 18AM 2 r+18Af 3 
V 1 (r) = 1 - — -oTT- , own9 . 



r 3 (Ar+3A/) 5 



(5) 



being A = 1/2(1— l)(l + 2). The source <S'(i, r) includes the derivative of the Dirac distribution, which appears in the 
process of building a single wave equation out of the seven even linearised wave equations: 



2(r - 2M)k 



S l = 



r 2 (A + l)(Ar + 3M) 



r(r-2M)_ /r /M 



r(A + 1) - 3M 3M[/°(r- - 2M) 2 



2U° 



r(Xr + 3M) 



6[r-r u (t)} 



(6) 



£7° = y/l - 2M/r u0 /(l - 2M/r u ) being the time com- different forms according to the initial conditions. For 
poncnt of the 4- velocity and K = Amy/ (21 + 1)tt. The the radial infall of a particle starting from rest at finite 
trajectory in the unperturbed SD metric r u (t) assumes distance r u0 , r u (t) is the inverse function of: 



t(r u ) 
2M 



= \ 1- 



2AI 



1-— (St) f^^ + Sarctanh 
r u0 \2MJ \2MJ 



/ 2M 2M\ 



\ 



1 



2M 



J 



1 - 



2M 



AM 
r u o 



\2M) 



3/2 



arctan 



- 1 



(7) 



r 



where E = y/l - 2M/r u0 . 

For an infalling mass from infinity at zero velocity, the 
energy radiated to infinity for all modes [37| is given by 
(in physical units): 



rad 



0.0104- 



M 



(8) 



while most of the energy is emitted below the frequency: 



fn 



0.08 



GM 



(9) 



Up to 94% of the energy is radiated between 8M and 
2M and 90% of it in the quadrupole mode. Obviously 
the above expressions would sensibly vary in case of an 
initial rclativistic velocity. 

For the SD black hole perturbations, there have been 
successive corrections of the basic equations, the last be- 
ing pointed out by Sago, Nakano and Sasaki (66[. Sago et 
al. have corrected the ZM equations (a minus sign miss- 
ing in all right-hand side terms) but introduced a wrong 
definition of the scalar product leading to errors in the 
coefficients of the energy- momentum tensor (67| . The ex- 
pressions herein correspond to those in (68| , where some 
of the errors of previously published literature on radial 
fall arc indicated. 



III. THE JUMP CONDITIONS 

It has been indicated by two different heuristic argu- 
ments (69l . FtqI ] that even metric perturbations, for radial 
fall in the RW gauge, should belong to the C° continuity 
class at the position of the particle. One argument [69( 
is based on the integration in r of the Hamiltonian con- 
straint, which is the tt component of the Einstein equa- 
tions (Eq. (C7a) in H3|); the other [7(| on the structure 
of selected even perturbation equations. Following the 
latter, Lousto and Nakano evince the continuity of the 
even perturbations and afterwards impose such continu- 
ity to derive the jump conditions on the wave function 
and its derivatives for I = 2, notably starting from the 
ZM equation. 

In this section, we instead provide an analysis vis a 
vis the jump conditions that the wave function and its 
derivatives have to satisfy for guaranteeing the continu- 
ity of the perturbations at the position of the particle. 
This approach [tT| is based on the solutions of the ZM 
equation, not on the equation itself. The jump conditions 
found herein are applicable to all modes. Incidentally, the 
jump conditions were previously mentioned by Sopuerta 
and Laguna [72| , Haas [73| , and Field et al. ukli in a dif- 
ferent context. In the forthcoming paper |64j |. the jump 
conditions for generic orbits are based on the analysis 
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of the odd and even wave equations and are determined 
despite the lack of continuity of the perturbations. 

In the radial case, only three independent perturba- 
tions survive, namely A'', H l = H 2 , H[. The inverse re- 



lations for those perturbations, as functions of the wave 
function and its derivatives, are given by (we drop hence- 
forth the I index): 



r ^ 6M 2 + 3MXr + A(A + l)r 2 / 2M\ K [/°(r-2M) 2 . 

A = — — )- J —y +1 r - — -i- otuF\~ ' ( 10 ) 



r 2 (Xr + 3M) \ r J ' (A + l)(Ar + 3M)r ' 

9M 3 + 9AM 2 r + 3 A 2 Mr 2 + A 2 (A + l)r 3 MI 2 - XMr + Ar 2 
H„ = : : v]> J \[> -L- 

r 2 (Ar + 3Af) 2 r(\r + 3M) 

„,^ T K [/°(r-2M)(A 2 r 2 + 2AA/r-3A/r + 3Af 2 ) c . K U°(r~2M) 2 r . . , 

( r - 2M ^-rr + r(A + l)(Ar + 3AA) 2 (A + l)(Ar + 3M) ' (U) 



_ Ar 2 - 3MAr - 3M 2 k U° z u (Ar + M) k U° z u r(r - 2M) 

1 (r - 2M) (Ar + 3M) '* '* r (A + 1) (Ar + 3Af ) (A+l)(Ar + 3M) ' [ ' 

where 8 = 5 [r — r u (t)] and 8' = 6' [r — r u (t)]. 



From the visual inspection of the ZM wave equation 
(|4]), it is evinced that the wave function ^ is of C -1 
continuity class, since the second derivative of the wave 
function is proportional to the first derivative of the Dirac 
distribution (incidentally, the concept of C _1 continuity 
class element, like a Hcaviside step distribution, may be 
pragmatically introduced as an element of a class of func- 
tions which after integration transforms into an element 
belonging to the C° class of functions). Thus, the wave 
function and its derivatives can be written as: 

# =9+(t,r) 6i + tf - (i,r) 2 , (13) 



of the particle, namely f(r)8'[r — r u {t)) = f\r u (t)S'[r — 
r u{t)\ — f'\r u (t)fi[ r ~ r u{t)}, has been used. 

The wave function and its derivatives are to be re- 
placed in Eqs. (|10M12[) . At the position of the particle, it 
is wished that the combination of all discontinuities ren- 
ders the perturbations of C° class. Thus, the coefficients 
of Oi must be equal to the coefficients of 02, while the 
coefficients of S and 8' must vanish separately. To this 
end, the perturbations are expressed in an implicit form: 

K= f 1 (r)V + f 2 (r)y tr + f 3 (r)5 , (18) 



$ r = #+0i +*-6 2 + (*+-*-) 8 , (14) 



"(15) 

*. t = *+ei + ^e 2 -(* + -*-)r u( 5 , (16) 



* >tr =*+ r e 1 +r^e 2 +(*+-^) *-(*+-*-) r f u 8' , 

" (17) 

where 9i = [r — r u (t)], and 02 = [r„(t) — r] are 
two Heaviside step distributions. For Eqs. (|15|17jl . a 
property of the Dirac delta distribution at the position 



# 2 = / 4 (r)tf + f 6 (r)% r + / 6 (r)*,rr + fr(r)8 + f s (r)5' , 

(19) 

Hi = f 9 (r)*, t + f 10 (r)*,tr + fn(r)8 + fa(r)5' , (20) 

where the definitions of the / functions may be easily 
drawn by visual inspection of Eqs. ()10ti 12[) . The jump 
conditions are obtained via Eq. (TT8"|) : 

(* + -^) ru =-f > ( 21 ) 
(*+-*") p =4^ - ( 22 ) 

viaEq. (Jigi: 
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(23) = _Ui^_?n + k(n^ 



(26) 



(Kr-Kr) r = T ( ~ * + k,r - M * 



and finally via Eq. (j2"0)) : 



k\ k k J d(^+ - . , T . . 

(24) (*+ - ^ = ^ - (*+ - *") r u 



=^_. (25) (/9-/io,)r u (^-^)-/n + ? (2?) 



The jump conditions set by Eqs. ([2l]). ([23]) and (J2SJ) 

are equivalent, as those set by Eqs. ([2"2"jl and (j2~4")l . The , , 

equivalence shows the consistency of the conditions and /^+ _ ^y- \ _ ^ 9 ~ ZAl (28) 

the compliance to the continuity requirement. The first ' r ru /io 

t, second r, second mixed t, r derivative jump conditions 

(coming from Hi and Hi) are given by: In explicit form, the jump conditions become: 



, . x kE\6M 2 + 3MXr u + X(X + l)r 2 u ] , 

(^+-^r _ ) = l - i-Hi (30) 

V r - r '^ (A + l)(2M-r„)(3M + Ar u ) 2 ' 1 J 

[3M 3 (5A - 3) + 6M 2 A(A - 3)r u + 3MA 2 (A - l)r 2 - 2A 2 (A + l)rj] 

1 ' rr < rrj ''" (A + l)(2A/-r u ) 2 (3M + Ar u ) 3 ' 1 J 



' a 

'»■« " (2M-r„)(3M + Ar„) ' 



(*j - * *l - - 7i»7 rcr," . n a , (32) 



_ kE 1 (3M 2 + 3M\r u - Ar 2 ) r M 

1 < tr ' tr) ^ (2M-r u ) 2 (3M + \r u ) 2 ' 1 ' 



Having established the jump conditions for any value 
of I at first order, we now turn to the presentation of the 
integration scheme. It is remarked that first derivatives 
of the jump conditions suffice to a scheme at first order. 



Each cell has an area of 2h 2 where \f~2h is the dimension 
of an edge of the cell. The evolution scheme starts with 
the initial data approach by Martel and Poisson [5(J . 



IV. THE NUMERICAL METHOD 

Our integration domain is discretised using a two- 
dimensional uniform mesh based on the coordinates t, 
r* as depicted in Fig. [TJ where -450 < r* jlM < 1000 
and < t < t en d, with tend greater than the falling time. 



A. The numerical scheme 

The jump conditions determined in the previous sec- 
tion are functions of the t, r variables, but for convenience 
in the following computations, they are transformed into 
functions of t, r* . The integration method considers cells 
belonging to two groups: for cells which are never crossed 
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FIG. 1: Numerical domain: a staggered 1+1-dimensional 
mesh in t,r* coordinates. The geodesic trajectory of the par- 
ticle is represented by the continued line 7. The shadowed 
cells are those crossed by the particle. 




FIG. 2: Cell crossed by the particle. The particle crosses the 
segment [a 6] at the point a and the segment [era] at the point 
b. 



(*+-*".)„ = [*,r.]« 



tb 



(36) 



(37) 



(38) 



(39) 



We maintain the superscript notation for which the 
wave function values between the world line and infinity 
are noted by the plus sign, and conversely by the minus 
sign for the values between the world line and the horizon 
(incidentally, such notation would also be applicable to 
the LPMP method). 

Using a first order expansion, we get a set of six nu- 
merical equations. For case 1, see Fig. they are: 



*-(t 6 -(h- e b ), r*) = - (h - efc) 



* lb 



(40) 



(41) 



*r = y-(t b -2h+e bl r* b ) = ^-{t a -h,r* b ) = *--h*;. 



= *+ (t a , r* a + t a ) = + e a *+ 



(42) 



(43) 



by the particle, the integrating approach is drawn by the 
LPMP method, whereas for cells which are crossed, we 
propose our strategy. 

There arc four sub-cases representing the different par- 
ticle trajectories inside the cell, see Figs. <J2] - [5]) . Wc 
define a, /?, 7, 5 the four vertices of the diamond cen- 
tred on <7 and initially consider the world line cross- 
ing the cell as displayed in Fig. ©. The trajectory 
crosses the line (3-5 at the point a and the line a-7 at the 
point b. We also define the shift e a and the lapse e b as 

e a = min jr* - r* 3 ,r* s - r*|, e b = min{i Q - t b ,t b - i 7 }, 

respectively. The jump conditions at the a, b coordinates 
express six analytical equations: 



*- = *-(t.,r:-(/i-e )) = *--(&-e )*-J , (44) 



I77 = *-{t a ,rl-2h+e a ) = *-(t ff ,r*-ft) = *T r ,| c 

(45) 

Our aim is the determination of the value of know- 
ing those of e a , e 6 , [*] a f) , [# r ] o6 and 
[* t ] t To this end, we subtract Eq. (gTJ) from Eq. (gDJ, 
Eq. (|4~4l from Eq. (1431) and obtain, respectively: 



= + [*] 6 + e 6 [* t ] 6 + ft 



* lb ' 



(46) 



(*+-*-)„ = [*]„ 



(34) 



(35) 



(47) 



Summing Eq. and Eq. (gBJ), Eq. and Eq. 

(|47p. and combining the results, we finally get: 
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*t=*e-% + n-Wa + l*} b -t a l*,r*} a + e b [*,t] b ■ 

, (48) 

We thus have obtained, without direct integration of 
the singular source, the value of the upper node. Further- 
more, the latter depends solely on analytic expressions, 
and obviously on the other values at the cell corners. 
Similar relations are found for the other three cases. For 
case 2 (Fig. [3]), we obtain: 

*i = *^-*^ + *? + [*]«-[*]6-e-[*,r.] + 66[*,t]i 5 

(49) 



for case 3 (Fig. 0J: 

*- = *--*"+*+- [*] a - e a [*. r .] a ; (50) 
for case 4 (Fig. |5j|: 

*£ = *fl-*} + *J + [*]a-M*,r.] a ■ (51) 

In a concise form, the four cases might be represented 
by a single expression (where ty a stands for or *5?~, 
and <3/ 7 stands for >]>+ or \J r ~): 



* a = *p-% + i!>}T{m a -B[*} b }-e a [*, r *} a + Be b [*, t } b , (52) 

where the upper sign holds when r* > r* and the lower when r* < r*, B = if the particle does not cross the a-j 
line, and B = 1 otherwise. 




FIG. 3: Cell crossed by the particle. The particle crosses the FIG. 4: Cell crossed by the particle. The particle crosses the 

segment [0-7] at the point b and the segment [a/3] at the point segment [a 8] at the point a. 

a. 



V. RADIATED WAVE FORMS AT INFINITY 
AND THE WAVE FUNCTION AT THE 
PARTICLE 

We confirm the existing wave forms previously pub- 
lished by Lousto and Price (49|, and Martel and Poisson 
(50| . In this section, the results from a code based on the 
LPMP method (full second order) are compared to those 
obtained by the -I- method (first order for the filled cells, 
as presented herein, and LPMP-like for empty cells). In 
spite of the partially different order of the two codes, the 
difference between the wave forms computed by the two 
methods is marginal. Incidentally, this occurs in absence 
of a recognised standard to which refer different numeri- 
cal approaches. 

Figs. |51 [7] show the I = 2 wave forms for the LPMP 



and -I- methods, for a particle falling radially, with zero 
initial velocity, from t u q/2M = 5 and r„o/2M = 20, re- 
spectively, in units u = t—r*, (2M/m)f . The initial data 
are given by the minimal (a = 1) initial data condition 
H2 = olK for both cases (the parameter a, introduced 
in [50j , measures the amount of radiation present on the 
initial hypersurface and for a = 1, the initial metric is 
conformally fiat). 

Figs. [SI show the logarithmic (absolute) difference 
for the I = 2 wave forms between the LPMP and -I- 
mcthods, for a particle falling radially, with zero initial 
velocity, from r„o/2M = 5 and r u o/2M = 20, respec- 
tively It is evinced that the normalised wave function 
amplitudes differ of an amount between 10~ 3 and 10~ 8 5 , 
in normalised units. 

We show also the behaviour of the wave function at the 
position of the particle, Figs. (fTUI [TT|) . The theoretical 



FIG. 5: Cell crossed by the particle. The particle crosses the 
segment \<J0\ at the point a. 
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FIG. 6: Radiated I — 2 wave form for a particle falling radially 
from r u o/2M = 5 with zero initial velocity, in units it = t—r*, 
(2M/m)4 r . The dashed and solid line represent the -I- and 
LPMP methods, respectively. 

jump condition, Eq. (|29j) . and the numerical result well 
overlap at the particle position, as it may be inferred from 
the logarithmic (absolute) difference plots, Figs. ([12113)1 . 
The discrepancy in normalised units is between 10 -3 and 
10" 6 - 5 . 



A. Discussion 
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FIG. 7: Radiated I = 2 wave form for a particle falling 
radially from r u a/2M — 20 with zero initial velocity, in units 
u = t — r*, (2M/m)$. The dashed and solid line represent 
the -I- and LPMP methods, respectively. 
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FIG. 8: Logarithmic (absolute) difference plot for a particle 
falling radially from r u o/2M = 5 with zero initial velocity, 
in units u = t - r* , (2M/m)f, between the -I- and LPMP 
methods (1 = 2). 



The features of the -T method can be summarised as 
follows: 

• For the grid cells crossed by the particle, direct and 
explicit integration of the wave equation, the source 
term with the associated singularities and the po- 
tential, is avoided. This determines also a faster 
code. 



• Reliability is improved, since analytic expressions 
totally replace the numerical expressions represent- 
ing the source. The terminology "source-free" is 
due to this feature. 

• A higher order scheme can be written, starting from 
a higher order Taylor expansion for Eqs. (|40H45p , 
and associated to this method tfjfj . 
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FIG. 9: Logarithmic (absolute) difference pfot for a particle 
falling radially from r u o/2M = 20 with zero initial velocity, 
in units u — t - r* , (2M/m)$, between the -I- and LPMP 
methods (1 = 2). 
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FIG. 10: Wave function (*+, dashed-dotted line, and 
dotted line) at the particle position for a particle falling ra- 
dially from r u o/2M = 5 with zero initial velocity (1 = 2), in 
units u = t — r* , (2M/m)^. The numerical jump condition 
[$]mim, solid line, and the theoretical jump condition [\t]tfc e , 
dashed line Eq. (|29p , are shown. 
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FIG. 11: Wave function (* + , dashed-dotted line, and 
dotted line) at the particle position for a particle falling ra- 
dially from r u o/2M = 20 with zero initial velocity (1 = 2), in 
units u = t— r* , (2M /m)^ . The numerical jump condition 
[fjnom, solid line, and the theoretical jump condition [^Jthe, 
dashed line Eq. (|29[) , are shown. 
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FIG. 12: Logarithmic (absolute) difference plot for a particle 
falling radially from r u o/2M = 5 with zero initial velocity, in 
units u — t — r* , (2M/m)f , at the particle position, between 
theoretical and numerical jump conditions (1 = 2). 



• The applicability of the method stretches out to 
generic orbits. It is our concern to apply the - 
I- method to circular and eccentric orbits, even if 
these orbits are not accompanied by perturbations 
of C° class. The rationale poses on the following 
consideration. Instead of using the continuity of 
the perturbations to get the jump conditions, we 



assume that the even and odd wave equations are 
satisfied by respectively R, being C^ 1 . Using 
this approach, we have obtained encouraging re- 
sults for an eccentric orbit f&Ij . 

• The disadvantage of this method is represented by 
the four possible trajectories that a particle may 
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FIG. 13: Logarithmic (absolute) difference plot for a particle 
falling radially from r u o /2M = 20 with zero initial velocity, in 
units u — t — r* , (2M/m)>t, at the particle position, between 
theoretical and numerical jump conditions (1 = 2). 

follow inside a given cell, instead of the three possi- 
ble paths of the LPMP method. A different labeling 
of the intersections between the particle world line 
and the cell, though, reduces the number of cases 
to three (65j . 

VI. CONCLUSIONS AND OUTLOOK 

We have presented a novel integration method in the 
time domain for the Zcrili-Moncrief wave equation at first 
order, for cells crossed by the particle world line. The for- 
ward time value in the upper node of the (t, r*) grid cell 
is obtained by the linear combination of the three pre- 



ceding node values and of the analytic jump conditions. 
Therefore, the numerical integration does not deal any 
longer with the source term, the associated singularities 
and the potential. The direct integration of the wave 
equation is circumvented. 

The wave forms at infinity confirm published results, 
and the wave function at the particle position shows a 
well-behaved pattern. Our final aim is the evaluation 
and the application of the back-action at the position of 
the particle, in a self-consistent manner. 

The indirect or source-free method has been developed 
to fourth order [65| and applied to generic orbits [gU . 

Beyond the scenario of EMRIs, other lines of inves- 
tigation may emerge. One concerns the formal relation 
between the al gori thms of the indirect and LPMP 
approaches |4!| l5C)| . Another considers applicability 
of the indirect method to all wave equations with a 
singular source term, and which the wave function of is 
constrained by a set of properties. 
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